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Abstract 

We  present  a  new  algorithm  based  on  Wiener-Hermite  functionals  combined  with  Fourier  collocation  to  solve 
the  advection  equation  with  stochastic  transport  velocity.  We  develop  different  stategies  of  representing  the 
stochastic  input,  and  demonstrate  that  this  approach  is  orders  of  magnitude  more  efficient  than  Monte  Carlo 
simulations  for  comparable  accuracy. 
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1  Introduction 


Most  of  the  effort  in  CFD  research  so  far  has  been  in  developing  efficient  algorithms  for  different  applications, 
assuming  an  ideal  input  with  precisely  defined  computational  domains.  With  the  field  reaching  now  some  degree  of 
maturity,  we  naturally  pose  the  more  general  question  of  how  to  model  uncertainty  and  stochastic  input,  and  how 
to  formulate  algorithms  in  order  for  the  simulation  output  to  reflect  accurately  the  propagation  of  uncertainty.  To 
this  end,  the  Monte  Carlo  approach  can  be  employed  but  it  is  computationally  expensive  and  it  is  only  used  as  the 
last  resort.  The  sensitivity  method  is  an  alternative  more  economical  approach,  based  on  moments  of  samples,  but 
it  is  less  robust  and  it  depends  strongly  on  the  modeling  assumptions  [1].  There  are  other  more  suitable  methods 
for  physical  applications,  and  there  has  already  been  good  progress  in  other  fields,  most  notably  in  seismology  and 
structural  mechanics.  A  number  of  papers  and  books  have  been  devoted  to  this  subject,  most  notably  with  the 
work  of  Ghanem  and  others,  e.g.  [2,  3,  4,  5,  6,  7,  8,  9].  On  the  theoretical  side,  the  work  of  Papanicolaou  and  his 
collaborators  has  addressed  many  aspects  of  wave  propagation  and  other  physical  processes  in  random  media,  see 
for  example  [10]. 

An  effective  approach  pioneered  by  Ghanem  &  Spanos  [5]  in  the  context  of  finite  elements  for  solid  mechanics  is 
based  on  a  spectral  representation  of  the  uncertainty.  This  allows  high-order  representation,  not  just  first-order  as 
in  most  perturbation-based  methods,  at  high  computational  efficiency.  It  is  based  on  the  original  ideas  of  Wiener 
(1938)  on  homogeneous  chaos  [11,  12]  and  the  Hermite  polynomials.  The  effectiveness  of  Hermite  expansions  was 
also  recognized  by  Ghorin  (1971)  [13]  who  employed  Wiener-Hermite  series  to  substantially  improve  both  accuracy 
and  computational  efficiency  of  Monte  Garlo  algorithms.  However,  the  limitation  of  prematurely  truncated  Wiener- 
Hermite  expansions  in  applications  of  turbulence  has  also  been  diagnosed  [14]. 

In  this  paper  we  consider  the  linear  advection  equation  with  a  random  transport  velocity  as  a  prototype  problem, 
and  as  a  first  attempt  to  model  uncertainty  in  propagation  phenomena  using  polynomial  chaos.  The  objective  is  to 
evaluate  the  Wiener-Hermite  polynomial  chaos  in  terms  of  its  accuracy  and  efficiency.  Given  the  simple  equation 
considered,  we  are  able  to  derive  an  exact  solution  for  the  mean  and  the  variance  of  the  stochastic  solution.  A 
new  element  in  this  paper  is  the  representation  of  the  stochastic  input  using  Karhunen-Loeve  expansions  both  for  a 
specified  covariance  kernel  as  well  as  a  kernel  constructed  explicitly  following  a  dynamical  systems  approach.  In  the 
simulations,  we  consider  both  Gaussian  and  log-normal  distributions  of  the  transport  velocity.  In  the  next  sections, 
we  first  formulate  the  polynomial  chaos  algorithm,  derive  the  exact  solutions,  and  subsequently  discuss  different 
ways  of  input  representation.  We  then  present  simulations  with  time-  and  space-dependent  transport  velocity,  and 
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we  conclude  with  a  brief  summary. 

2  Polynomial  Chaos:  Wiener-Hermite  Expansions 

We  consider  the  stochastic  advection  equation 

^  +  F(x,  0)  ^  =  0  V(t,  x)  G  [0,  T]  X  I?  =  [-1, 1] 

^  u{t  =  Q,x)  =  g{x)  Vx  G  I? 

Periodic  boundary  conditions 

Here  0  is  a  random  parameter.  We  will  examine  cases  with  both  time-independent  as  well  as  time-dependent  transport 
velocity. 

We  first  assume  that  the  transport  velocity  V {x,  &)  is  only  space-dependent  and  is  a  given  random  process  with 
mean  value  and  variance  satisfying 

V{x)  =  E{V{x,9)}  =  1  and  E{\V{x,9)\^}  <  oo  (2) 

where  the  covariance  function 


C{x,y)=E{V{x,9)V{y,9)} 


(3) 


is  continuous.  Its  orthogonal  decomposition  is  obtained  by  the  Karhunen-Loeve  expansion 

M  M 


V(x,9)  =  y(x)  ^  \f^fk{x)^k{9)  =  y^gfc(x)^fc(6>), 


(4) 


k=l  k=0 

where  are  Gaussian  variables,  and  M  is  the  dimension.  The  eigenvalues  Afc  and  corresponding  eigenfunctions  fk 
satisfy  the  Helmholtz  integral  equation 


C{x,y)fk{y)dy  =  \kfk{x). 


(5) 


iv 


By  assuming  that  u(t,  x,  9)  is  a  second-order  process,  that  is  u(t,  x,  9)  has  a  finite  variance,  the  theorem  of  Cameron 
&  Martin  [15]  guarantees  that  the  solution  u{t,x,9)  can  be  represented  in  terms  of  polynomials,  i.e. 


i{t,  X,  9)  —  agPo  +  +  E  E 

=  l  ii—1  i2—l 

OO  il  l2 

+  E  E  E  +  •  ■  • 


(6) 


*1  =  1  *2  =  1  *3  =  1 

In  the  above  equation,  rn(^ii , C*2 j  •  ■  •  Ain)  denotes  the  polynomial  chaos  of  order  n  in  the  Gaussian  variables 
(C*l  1  C*2  1  ■  5  An  )  ■ 
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The  general  expression  for  generating  these  polynomials  is  given  by 


i9”e 


(7) 


where  the  vector  ^  consisting  of  the  n  random  variables  ,  ^i2j  •  ■  •  >  '?i„  )• 

Having  equation  (7),  the  decomposition  in  equation  (6)  can  be  thought  of  as 

Gaussian  part 

/ - ^ 

OO 

u{t,x,9)  =  aoTo  +  ^  aiiri(^i^(6»)) 

OO  00  21  22 

" - V - " 

non-Gaussian  part 

For  a  Gaussian  process,  only  the  first  term  is  present  and  thus  this  expansion  reduces  to  the  Karhunen-Loeve 
representation  of  the  same  process  [16]. 

Equation  (6)  can  be  re-written  in  the  more  familiar  form 

p 

u{t,  x,e)  =  '^  a,{t,  a;)Ti(|(6l)),  (8) 

i=0 

where  there  exists  one-to-one  correspondence  between  the  and  '£'^(^(0))  as  well  as  between  the 

coefficients  ai{t,x)  and  The  total  number  of  polynomial  chaoses  {P  +  1)  used  in  the  expansion  is 


(M  +py. 
M\p\ 


where  M  and  p  refers  to  the  dimension  of  the  stochastic  input  and  the  order  of  the  polynomial  chaos,  respectively. 

In  equation  (8),  ai{t,x)  are  the  deterministic  coefficients  and  is  a  trial  basis  in  the  space  of  random 

variables.  This  basis  is  the  set  of  generalized  multi-dimensional  Hermite  polynomials  in  the  quantity  ^{0).  This 
representation  can  be  refined  along  the  random  dimension  either  by  adding  more  random  variables  (dimension  M) 
or  by  increasing  the  maximum  order  of  polynomials  included  in  the  polynomial  chaos  expansion  (order  p) . 

If  the  solution  is  known,  then  the  coefficients  ai  in  the  expansion  (8)  can  be  obtained  from 

<  u{t,x,9),'^i{i{0))  > 


ai{t,x)  = 


<vI/,(e(0)),T,(e(0))  > 


(9) 


where 


Jn 


(10) 


is  the  inner  product  (statistical  average). 
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Substituting  the  above  in  the  advection  equation  (1),  it  yields 

Q  -P  M  ^  P 

—  a,{t,  x)'^i{^{e))  +  (Y  Y  a;)«'i(^(6»))  =  0.  (11) 

i— 0  j—0  i—0 

Assuming  that  we  discretize  in  space  with  a  spectral  method  we  write 

N 

a^{x,t)  =Y'^ik{t)^k{x)  (12) 

k=0 

where  <I)fe(x)  denotes  the  trial  basis  in  space;  here,  we  employ  the  Fourier  exponentials,  i.e.  d)fc(x)  =  .  Upon 

substitution,  we  obtain 

P  N  f)  (j.\  p  M  N  a(f,  (  \ 

i—O  k—0  i=0  j—0  k—0 

We  now  apply  a  standard  Galerkin  projection,  first  in  space  x 

EE  j  ^k{x)(^i{x)dx+ 

i=0  k=o  ''P 

P  M  N  „  /  X 


EEE^l(^)^*(^(^))'^*'=(^)  /  <^i{x)dx  =  0 

„• _ r»  _ n  j _ n  'J  T) 


2=0  i=0  k—0 


for  I  =  0,1,2,...  ,1V 


We  denote  by 


the  mass  matrix,  and  by 


Iki  =  /  ^k{x)^i{x)dx 
Jv 


Ijki  =  J  gj{x) — — 9i{x)dx, 


the  derivative  matrix,  to  arrive  at 


EE'«^**K>))+EEE  Ijkiaik{t)^j{6)'i>i{^{0))  =  0 

2=0  k—0  i=0  j—0  k—0 

for  1  =  0, 1,2,...  ,1V 


P  M  N 


By  analogy  with  the  deterministic  approach,  we  use  a  Galerkin  projection  for  the  random  variable  as  well,  i.e.  we 
multiply  the  system  (15)  by  'I'm(C(^))  and  take  the  statistical  average  defined  in  equation  (10).  The  orthogonality 
property  of  '£'^(^(0))  is  employed  to  arrive  at 


Y^ki 


P  M  N 

^jkl  Q^2fc(^)  ^ —  0 

2=0  ^=0  k—0 


for  /  =  0, 1,  2, . . .  ,  and  m  =  0, 1,  2, . . .  ,  P 
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Once  the  coefficients  aik  have  been  evaluated,  the  stochastic  solution  can  be  obtained  from 

p  N 

u{t,x,9)  =  (17) 

z^O  k—0 

where  the  superscript  n  here  denotes  time  discretization. 

We  close  the  section  by  pointing  out  that  the  previous  analysis  remains  unchanged  in  the  case  of  V (x,  9)  being 
a  non-Gaussian  process.  In  that  case,  the  polynomial  chaos  expansion  is  used  to  represent  the  process  instead  of 
the  Karhunen-Loeve  expansion.  Consequently,  the  term  <  j'i> i'i> m  >  replaces  <  >  in  equation  (16).  A 

more  general  approach  that  involves  different  trial  bases  from  the  Askey  scheme  of  orthogonal  polynomials  has  been 
developed  in  [17]. 

3  Exact  Stochastic  Solutions 


If  the  transport  velocity  V  in  equation  (1)  is  a  random  function  of  the  time  t  alone,  the  initial  value  problem  can  be 
solved  exactly  by  the  method  of  characteristics  or  by  a  change  of  variables.  We  demonstrate  the  latter  method  here. 
We  remove  the  random  positive  variable  V  (t)  from  the  equation  by  a  change  of  variable  t  to  r  as  follows 

dr  =  V  (t)  dt 
T  =  0  when  t  =  0. 


The  advection  equation  (1)  is  then  reduced  to 


du  du 
dr  dx 

M(a;,0)  =  g{x). 


The  solution  to  the  above  equation  is  simply  given  by 


(18) 


(19) 


u{x,t)  =  g{x  -  r). 


(20) 


while  the  solution  for  t  is  obtained  by  solving  equation  (18) 


^{t)  =  /  V (s)ds  =  At 


Vi. 


i=i 

In  order  to  evaluate  the  sum  of  the  set  of  random  variables  Vj  in  (21),  we  assume  Vj  to  be 

vi  =  a^i 

Uj  =  cvi-i  +  af^i  for  i  =  2,3, . . .  ,Q 


(21) 


(22) 
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where  a,  c  determine  the  degree  of  correlation  of  the  input,  and  /^  =  1  —  to  ensure  that  each  Vi  has  the 
variance  a.  Also,  ^i,  ^2)  ■  •  ■  j  Cq  ^re  a  set  of  independent  normal  random  variables  of  variance  unity,  i.e., 

1  -i! 

?i  =  ^e  2  Vz=l,2,...,g 

V  27r 

We  obtain  the  following  three  cases  : 

Case  I  -  Fully  Correlated  Input:  For  c  =  1  and  /  =  0  the  input  is  fully  correlated,  and 

Q 

Vj  =  aQ^i  =  aQ^ 
i=i 

Case  II  -  Mutually  Independent:  For  c  =  0  and  /  =  1  the  input  is  mutually  independent,  and 

Q 

i=i 


Case  III-  Partially  Correlated  Input::  For  general  values  of  c,  we  obtain 


i=i 


Letting  t  =  QAt,  we  obtain 


Q 


r{t)  =  At  '^3  = 

i=i 


at^  +  Vt 

_  _  aV tAt^  +  Vt 

i[2At  -  2A^{1  -  +  Vt  as  (c  ^  1,  At  ^  0) 

We  can  now  write  the  solution  to  equation(20)  as 

u{x,  t)  =  u{x,  r(t))  =  g{x  -  T{t)) 

and  obtain  its  mean  by  taking  the  expectation  value  of  equation  (25)  with  the  normal  distribution  (23): 

(x  —  xo  —  VtY 


1  f°° 

u{x,t)  =  —^ —  /  9{xo)e 


-\/^i 


TTCLCT  J  —  oo 


2a2(j^  dxQ  =  sin7r(x+  1  —  Vt)e 


where  cr  is  defined  from 


=  < 


\  t^  for  fully  correlated  case 

{At)t  for  mutually  independent  case 

2A\t  —  A{1  —  for  partially  correlated  case 


same 


(23) 


(24) 


(25) 


(26) 


(27) 
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The  variance  of  the  random  variable  u{x^  t)  is  obtained  from 


Var{u)  =  E{v?}  —  v? , 

which  together  with  the  normal  distribution  (23)  yields 

Var{u{x,t))  =  ^[1  -  cos2^(cc  -  =  i(i  _  +  cos27r(a;  +  1  - 

where  g{xQ)  was  set  to  sin7r(a:o  +  1)  in  order  to  obtain  equation  (28). 


(28) 


4  Representation  of  Stochastic  Inputs 


We  consider  different  representations  of  the  stochastic  inputs  corresponding  to  covariance  kernels  specified  or  con¬ 
structed  explicitly  following  a  dynamical  systems  approach. 

4.1  Specified  Covariance  Kernels 

In  the  case  of  time-dependent  but  space-independent  transport  velocity,  i.e. 

V{t,x,ijj)  =  V{t,uj) 


we  use  the  exponential  kernel 


C{t,s)  =  exp{--^^-^}, 


where  6  is  a  parameter  which  has  the  same  dimension  as  time  t  and  is  termed  correlation  length.  It  expresses  the 
rate  at  which  the  correlation  decays  between  two  time  instants  of  the  process.  For  the  kernel  under  consideration, 
the  analytical  solution  to  the  integral  equation  (5)  can  be  found  in  the  following  fashion: 

Let  T  =  [0,T]  where  T  designate  the  final  time,  and  set  c  =  ■^.  The  integral  equation  (5)  becomes 


[  e  ^1*  ^'/fe(s)ds  =  Xkfk{t), 
Jo 


(29) 


which  can  be  written  as 

f  e-‘^^*-^^fk{s)ds  +  r  fk{s)ds  =  Xkfk{t), 

Jo  Jt 

Since  the  integrands  in  (30)  are  continuous,  differentiating  (30)  twice  with  respect  to  t  yields 


(30) 


Xf\t)  =  {-2c  +  c^X)f{t), 


(31) 


where  the  subscript  k  was  omitted.  The  boundary  conditions  associated  with  (31)  are 


r  /(o)-c/(o)  =  o 

1  /(r)  +  c/(r)  =  o 

Making  the  substitution 

9  2c  —  c^A 


the  boundary  value  problem 

r  /(t)+u;VW  =  0 

<  /(0)-c/(0)  =  0 

[  /(r)  +  c/(T)  =  o 


(32) 


(33) 


(34) 


follows.  Solutions  to  (34)  are  given  by 


f{t)  =  acos{uit)  +  Psin{u!t), 


(35) 


with  dispersion  relation 


(w^  —  c^)  tan  wT  —  2cw  =  0, 


from  which  lu  can  be  determined.  The  normalized  eigenfunctions  are 


fit) 


:Cos{ujt)  +  sin{ijjt) 


1 ..  ^,sin{2ujT) 


-(i  +  -)r+(--i): 


4w 


—  (1  -  cos(2wT)) 


(36) 


(37) 


When  the  transport  velocity  is  space-dependent, 


V {t,  =  V {x,  w) 


we  use  a  similar  kernel  of  the  form  employed  by  Ghanem  &  Spanos  [18] 


C{x,y)  =  exp{-  ^  } 

A  similar  analysis  as  before  produces  analytical  forms  for  the  eigenvalues  and  eigenfunctions. 

4.2  Dynamical  Systems  Approach 

We  also  consider  different  ways  for  correlations  to  represent  stochastic  inputs,  which  can  be  useful  to  represent 
discrete  input.  Specifically,  we  consider  the  following  two  processes: 

Uk  =  cuk-i  +  af^k,  (38) 
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(39) 


Uk  =  2  (^fc+i  +  Uk-i)  +  afik- 

The  first  process  (38)  is  autoregressive  of  order  1  and  corresponds  to  the  Markov  process.  The  second  process  is 
termed  bilateral  autoregressive  process  due  to  its  backward  and  forward  dependence.  The  constant  c  is  assumed  to 
be  0  <  c  <  1  in  order  to  ensure  that  the  processes  is  of  finite  variance.  Here,  ^k  is  a  random  variable  of  mean  zero 
and  variance  one,  and  /  is  a  constant  to  be  determined  such  that  for  the  given  values  of  a  and  c  the  variance  of  the 
process  is  equal  to  a^. 

We  can  construct  numerically  the  covariance  kernel  and  subsequently  extract  the  eigenvalues  and  eigenfunctions 
needed  for  the  input.  Alternatively,  we  can  write 

N-k+l  N 

Uk  =  ^  Oii^k+t-1  +  ^  Oit^N-k+1  (40) 

i—1  i—N  —  k-\-2 

and  solve  for  A  simple  calculation  based  on  projection  and  orthogonality  of  the  Gaussian  random  variable 
yields  the  following  expression  for  the  mean  value 

\i-i\ 

E^Ui'Uj'\  ^  ^  ^k0^k+\i—j\  T  ^  ^  (|z— j|— fc)  ■  (41) 

fc=l  fe=l 

In  figures  1  and  2  we  give  the  eigenfunctions  and  compare  the  eigenspectra  between  a  unilateral  and  a  bilateral 
autoregressive  process.  Periodic  boundary  conditions  were  assumed  for  these  calculations. 


Figure  1:  Eigenfunctions  and  eigenspectra  corresponding  to  c  =  0.5  and  a  =  1. 
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Figure  2:  Comparison  between  unilateral  and  bilateral  eigenspectra. 

5  Numerical  Results 

We  consider  the  stochastic  advection  equation  with  initial  condition 

5(a:)  =  sin7r(l  +  a:)  xe[— 1,1]. 

In  the  following  we  use  Fourier  collocation  to  treat  the  deterministic  contributions.  We  also  perform  Monte  Carlo 
simulations  to  compare  with  the  polynomial  chaos  results.  For  the  latter  we  found  it  necessary  to  use  a  Lagrangian 
approach  in  solving  the  stochastic  advection  equation.  That  is,  we  have  implemented  a  solution  based  on  charac¬ 
teristics  that  avoids  global  spectral  interpolation.  In  contrast,  following  an  Eulerian  approach  leads  to  erroneous 
results  due  to  the  lack  of  smoothness  in  the  stochastic  solution.  Specifically,  the  Lagrangian  and  Eulerian  numerical 
solutions  are  initially  very  close  to  each  other  but  for  longer  time  integration  errors  in  the  Eulerian  approach  render 
the  solution  erroneous. 

In  the  following,  we  consider  several  cases  corresponding  to  different  types  of  stochastic  input  representation.  We 
employ  both  ad  hoc  type  covariance  kernels  as  well  as  kernels  constructed  based  on  a  dynamical  systems  approach. 


11 


5.1  Space-Independent  Transport  Velocity 

First,  we  assume  that  V{t,x,9)  =  V{t,9)  and  employ  the  kernel 

C{s,t)  =  exp{--^^^} 

for  the  covariance  matrix.  We  also  perform  a  companion  Monte  Carlo  simulation;  this  kernel  corresponds  to  a 
unilateral  autoregressive  process,  as  defined  previously.  For  comparisons,  we  also  obtain  the  analytical  mean  solution 
and  variance. 

For  the  fully  correlated  case,  figures  3  and  4  show  the  results  for  the  mean  and  variance,  respectively.  Good 
agreement  between  the  solution  obtained  from  the  polynomial  chaos  simulation,  the  Monte  Carlo  simulation,  and 
the  exact  formulas  is  obtained. 


Figure  3:  Mean  solutions  for  the  fully  correlated  case,  after  t  =  1  ,  t  =  3  and  t  =  5.  The  mean  solution  decays 
exponentially  with  time. 

For  the  partially  correlated  case  and  the  mutually  independent  case  similar  conclusions  to  the  fully  cor¬ 
related  case  can  be  drawn.  Indeed,  as  shown  by  the  figures  5  and  6,  also  figures  7  and  8,  good  agreement  between 
polynomial  chaos  simulation,  Monte  Carlo  and  exact  formulas  is  established  for  both  mean  and  variance.  It  is  worth 
mentioning  that  in  the  case  of  partial  correlation,  up  to  M  =  6  stochastic  dimensions  were  required  for  the  polynomial 
chaos  simulation  in  order  to  agree  with  the  exact  formula. 
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Figure  4:  Variance  of  the  solution  for  the  fully  correlated  case,  after  t  =  3. 


Figure  5:  Mean  solution  for  the  partially  correlated  case,  c  =  0.5,  after  t  =  5;  M  =  6  and  p  =  4. 
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Figure  6:  Variance  of  the  solution  for  the  partially  correlated  case,  c  =  0.5,  after  t  =  5;  M  =  6  and  p  =  4. 


Figure  7:  Mean  solution  for  the  mutually  independent  case,  c  approaches  0,  after  t  =  3;  M  =  2  and  p  =  4. 
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Figure  8:  Variance  of  the  solution  for  the  mutually  independent  case,  c  approaches  0,  after  t  =  3;  M  =  2  and  p  =  4. 

5.2  Time-Independent  Transport  Velocity 

We  replace  now  the  transport  velocity  V{t,x,6)  by  V{x,9)  and  assume  that  the  stochastic  transport  velocity  is 
represented  by  the  covariance  matrix 

C{x,y)  =  exp{--^^-^}. 

An  algorithm  similar  to  the  one  used  for  time-dependent  transport  velocity  was  employed.  For  a  length  scale  of 
A  =  1,  for  the  fully-correlated  case  after  a  time  of  t  =  1  we  observed  that  the  mean  solution  diverges.  The  cause  is 
attributed  to  the  non-periodicity  of  the  specified  exponential  kernel.  To  circumvent  this  difficulty,  we  constructed  a 
periodic  covariance  kernel,  generated  numerically  as  discussed  in  section  3;  see  figure  9. 

We  obtain  the  results  presented  in  figures  10  and  11,  after  integrating  to  time  t  =  10.  Here  again  agreement 
between  polynomial  chaos  and  Monte  Carlo  simulation  of  50,000  realizations  is  established. 

5.3  Non-Gaussian  Input 

Lastly,  we  consider  a  case  where  the  transport  velocity  is  not  Gaussian,  but  is  given  by  the  log-normal  distribution 

L 

V(x,0)  =  e^  h=h  +  Y,hi^i 
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Figure  10:  Time-independent  case:  Mean  velocity  after  t  =  10. 
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both  methods  as  there  is  no  exact  solution  available  for  this  case.  The  accuracy  of  the  approximation  in  the  mean 
is  higher  than  in  the  variance  but  overall  similar  conclusions  are  valid  for  the  lognormal  distribution  as  well.  In 
[17]  a  more  systematic  study  for  several  different  distribution  functions  is  performed,  and  appropriate  orthogonal 
polynomial  functionals  are  introduiced  as  trial  basis. 
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LcgmvnnQidi f4  bulton  hi  -  Z  p-  1  rilti  +  o^JfZ  -  -0.<^  ^  ~  ^ 


Figure  12:  Lognormal  distribution  mean  solution  corresponding  to  hi  =  /12  and  h  =  —0.5  after  t  =  1.5 


Lognormal  distribution  kl=2  p=4  with  oijj=-(a^  +  a^)/2  =  -0.5 


Figure  13:  Lognormal  distribution  variance  corresponding  to  hi  =  ft,2  and  h  =  —0.5  after  t  =  1.5 
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6  Summary 


We  have  developed  a  stochastic  spectral  method  to  solve  the  advection  equation  with  a  transport  velocity  described 
by  a  Gaussian  or  a  log-normal  distribution.  It  employs  Wiener-Hermite  functionals,  which  are  Hermite  polynomials  of 
a  Gaussian  random  variable  which  in  turn  depends  on  a  random  parameter  0.  More  general  distributions  can  also 
be  handled  with  Wiener-Hermite  expansions  but  we  have  found  in  ongoing  work  [17]  that  better  representations  with 
other  trial  bases  from  the  Askey  scheme  of  polynomials  lead  to  faster  convergence  compared  to  Hermite  expansions. 

The  Wiener-Hermite  expansion  exhibits  exponential  convergence,  which  we  verified  for  the  space-independent  case 
for  which  an  exact  solution  was  obtained.  To  realize  exponential  convergence  the  stochastic  input  has  to  be  accurately 
represented  but  this  depends  on  the  correlation  length  of  the  stochastic  input.  In  the  case  of  poorly  correlated  input 
a  very  large  number  of  Karhunen-Loeve  terms  are  required  to  represent  the  input  and  fast  convergence  is  not 
achievable.  The  relatively  poor  resolution  properties  of  Hermite  expansions,  compared  to  other  spectral  polynomials, 
are  well  documented  in  the  literature  [19,  20].  However,  re-scaling  procedures,  as  done  in  [21],  can  be  applied  or  a 
change  of  the  trial  basis  from  the  Askey  scheme  can  be  employed,  as  demonstrated  in  [17],  to  accelerate  convergence. 
Unfortunately,  the  exact  resolution  requirements  are  problem-dependent  and  there  is  not  enough  experience  at  the 
moment  in  order  to  come  up  with  practical  rules. 

As  regards  efficiency,  we  have  compared  with  corresponding  Monte  Garlo  simulations.  Again,  for  the  exact 
solution  the  standard  Monte  Garlo  approach  we  employed  requires  200,000  realizations  to  achieve  errors  in  the  mean 
solution  of  order  of  10“^.  The  same  accuracy  can  be  achieved  with  a  10-term  polynomial  expansion  leading  to  a 
speed-up  factor  of  20,000.  Glearly,  better  Monte  Garlo  approaches  with  accelerated  convergence  may  reduce  that 
speed-up  but  it  will  still  be  more  than  a  factor  of  1,000.  However,  as  we  discussed  above,  for  stochastic  input  that 
is  very  weakly  correlated  the  polynomial  chaos  approach  requires  a  Karhunen-Loeve  expansion  with  many  terms. 
This  may  increase  significantly  the  dimensionality  and  thus  the  computational  complexity  of  the  polynomial  chaos 
approach. 
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